#ifndef BIRTNORM_HPP
#define BIRTNORM_HPP

#include "rtnomNC.hpp"
#define PHI13 0.4307
#define Lambda 0.68
class Birtnorm
{
	public:
		Birtnorm(void)
		{
			rt= new rtnormNC();
			boost::random::normal_distribution<> No;
			dn= new boost::math::normal_distribution<>(0,1); 
			rnorm= new RandomG::Random<boost::random::normal_distribution<> >(No);
			boost::random::uniform_01<> U;
			runif=new RandomG::Random<boost::random::uniform_01<> >(U);
		}
		~Birtnorm()
		{
		}
		
		mat operator()(mat Sig, mat a)
		{
			double rho=Sig(0,1)/(sqrt(Sig(0,0)*Sig(1,1)));
			double nu2=1-rho*rho;
			double a_1=a(0,0)/sqrt(Sig(0,0));
			double a_2=a(1,0)/sqrt(Sig(1,1));
			
			bool acc;
			bool exch=false;
			mat X(2,1);
			if(a_1<a_2)
			{
				double temp=a_2;
				a_2=a_1;
				a_1=temp;
				exch=true;
				temp=Sig(0,0);
				Sig(0,0)=Sig(1,1);
				Sig(1,1)=temp;
			}
			if(((rho>=0) & ((rho*a_1-a_2)>=0)) | ((rho<0) & (a_1< -PHI13)))//S+
			{
				acc=true;
				while(acc)
				{
					X(0,0)=(*rt)(0,a_1,pINF,1); 
					X(1,0)=sqrt(nu2)*((*rnorm)())+rho*X(0,0);	
					if(X(1,0)>a_2)
					{
						acc=false;
					}
				}

									
			}else if((rho < 0) & (rho*a_1-a_2<=0) & a_1> -PHI13){//S-
				
				acc=true;
				while(acc)
				{
					X(0,0)=(*rt)(rho*a_2,a_1,pINF,sqrt(nu2)); 
					double r1=-(rho*X(0,0)-a_2)/sqrt(nu2);
					double r2=(rho*a_1-a_2)/sqrt(nu2);
					double prob=(cdf(*dn,-r1)/pdf(*dn,r1))/c(r2);
					double u=(*runif)();
					if(u<prob)
					{
						acc=false;
						X(1,0)=(*rt)(rho*X(0,0),a_2,pINF,sqrt(nu2));
					}
				}

			}else if((rho<0) & (rho*a_1-a_2)>0 &(a_1>PHI13)){//M+
				acc=true;
				double o1=cdf(*dn,a_2/rho)-cdf(*dn,a_1);
				double o2=0.5*sqrt(nu2)*exp(-0.5*a_2*a_2)*cdf(*dn,-a_2*sqrt(nu2)/rho);
				double prob=o1/(o1+o2);
				double u1=(*runif)();
				if(u1<prob)
				{
					while(acc)
					{
						X(0,0)=(*rt)(0,a_1,a_2/rho,1);
						X(1,0)=sqrt(nu2)*((*rnorm)())+rho*X(0,0);	
						if(X(1,0)>a_2)
						{
							acc=false;
						}
					}
				}else{
					while(acc)
					{
						X(0,0)=(*rt)(rho*a_2,a_2/rho,pINF,sqrt(nu2));
						double u=(*runif)();
						double r1=-(rho*X(0,0)-a_2)/sqrt(nu2);
						double prob1=sqrt(2/PI)*(cdf(*dn,-r1)/pdf(*dn,r1));
						if(u<prob1)
						{
							acc=false;
							X(1,0)=(*rt)(rho*X(0,0),a_2,pINF,sqrt(nu2));
						}
					}
				}
				
			}else{//M-
				                                                                 
				acc=true;

				double theta=rho*(a_2+Lambda*sqrt(nu2));
				double tau1=cdf(*dn,-a_2/rho);
				double tau2=sqrt(nu2)/sqrt(2*PI)*(cdf(*dn,(a_2/rho-theta)/sqrt(nu2))-cdf(*dn,(a_1-theta)/sqrt(nu2)));
				tau2*=exp((theta*theta-a_2*a_2-2*Lambda*sqrt(nu2)*a_2)/(2*nu2))*d((rho*a_1-a_2)/sqrt(nu2));
				double prob=tau1/(tau1+tau2);
				double u1=(*runif)();
				if(u1<prob)
				{
					while(acc)
					{
						X(0,0)=(*rt)(0,a_2/rho,pINF,1);
						X(1,0)=sqrt(nu2)*((*rnorm)())+rho*X(0,0);
						if(X(1,0)>a_2)
						{
							acc=false;
						}	
					}
				}else{
					while(acc)
					{
						X(0,0)=(*rt)(theta,a_1,a_2/rho,sqrt(nu2));
						double u=(*runif)();
						double prob1=exp(Lambda*(a_2-rho*X(0,0)/sqrt(nu2)))*cdf(*dn,-(a_2-rho*X(0,0))/sqrt(nu2))/pdf(*dn,(a_2-rho*X(0,0))/sqrt(nu2))/d((a_2-rho*a_1)/sqrt(nu2));
						if(u<prob1)
						{
							acc=false;
							X(1,0)=(*rt)(rho*X(1,0),a_2,pINF,sqrt(nu2));
						}
					}
				}
			}

			X(0,0)*=sqrt(Sig(0,0));
			X(1,0)*=sqrt(Sig(1,1));
			if(exch)
			{
				double temp=X(0,0);
				X(0,0)=X(1,0);
				X(1,0)=temp;

			}
			return X;			

		}
		double c(double x)
		{
			double y;
			if(x>=0)
			{
				y=sqrt(PI/2);
			}else if(-(double)1/x<sqrt(PI/2)){
				y=-(double)1/x;
			}else{
				y=sqrt(PI/2);
			}
			return y;
		}
		double d(double x)
		{
			double y;
			double chi=exp(-Lambda*x)*cdf(*dn,x)/pdf(*dn,-x);
			if(chi<sqrt(PI/2))
			{
				y=chi;
			}else{
				y=sqrt(PI/2);
			}
			return y;
		}
	private:
		rtnormNC *rt;
		RandomG::Random<boost::random::normal_distribution<> > *rnorm;
		boost::math::normal_distribution<> *dn;
		RandomG::Random<boost::random::uniform_01<> > *runif;
};
#endif
